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ABSTRACT 

We consider the physical conditions under which supermassive black holes could have formed inside 
the first galaxies. Our SPH simulations indicate that metal-free galaxies with a virial temperature ~ 10'' 
K and with suppressed H2 formation (due to an intergalactic UV background) tend to form a binary 
black hole system which contains a substantial fraction (> 10%) of the total baryonic mass of the host 
galaxy. Fragmentation into stars is suppressed without substantial H2 cooling. Our simulations follow 
the condensation of 5 x lO^M© around the two centers of the binary down to a scale of ^ 0.1 pc. 
Low-spin galaxies form a single black hole instead. These early black holes lead to quasar activity before 
the epoch of reionization. Primordial black hole binaries lead to gravitational radiation emission at 
redshifts z ^ 10 that would be detectable by LISA. 

Subject headings: black hole physics — cosmology: theory — galaxies: formation — hydrodynamics — 
quasars: general 



L INTRODUCTION 

Supermassive black holes (SMBHs) are believed to pro- 
vide the power source of quasars via the accretion of sur- 
rounding gas (Salpeter 1964; Zeldovich 1964; Lynden-Bell 
1969; Rees 1984). They were able to form already at early 
cosmic times, as implied by the recent discovery of quasars 
at redshifs z > 6 (Becker et al. 2001; Djorgovski et al. 
2001; Fan et al. 2002). The existence of SMBHs with 
inferred masses of ^ IO^Mq, less than a billion years af- 
ter the big bang, provides important constraints on any 
formation scenario (Haiman & Loeb 2001). 

The gas physics involved in the formation of SMBHs is 
still not well understood (Rees 1984, 2002; Loeb & Rasio 
1994; Barkana & Loeb 2001). In this work, we attempt 
to simulate this process under the simplest and best pre- 
scribed set of initial conditions; those dictated by the early 
universe. The earliest galaxies are simple in that they are 
the first condensations of gas to grow out of the seed inho- 
mogeneities in the early universe. The composition of the 
primordial gas is determined by big bang nucleosynthesis 
(Buries, NoUett, & Turner 2001, and references therein) 
and any primordial magnetic fields are not expected to be 
dynamically significant. 

Previous numerical simulations of early luminous struc- 
tures (Bromm, Coppi, & Larson 1999, 2002; Abel, Bryan, 
& Norman 2000, 2002; Nakamura & Umemura 2001) have 
focused on the formation of stars inside the very first 
gaseous objects with a mass, ~ lO^M©, just above the cos- 
mological Jeans mass. Due to their low virial temperature 
(hundreds of K), fragmentation into stars is possible inside 
these objects through the formation of molecular hydro- 
gen, H2, which cools efficiently via rotational-vibrational 
transitions even at these low temperatures. However, a 
relatively modest UV flux is sufficient to photo-dissociate 
the fragile H2 molecules and suppress their role in cooling 



the gas (Haiman, Rees, & Loeb 1997; Ciardi et al. 2000; 
Haiman, Abel, & Rees 2000). A destructive flux of UV 
photons could be produced by a small, early population of 
stars. 

Here we consider dwarf galaxies with virial temperatures 
^ lO'' K inside of which atomic cooling is effective (Oh & 
Haiman 2002). These galaxies comprise a substantial pop- 
ulation at z ~ 10 (Barkana & Loeb 2001) that is expected 
to dominate the reionization of hydrogen at z > 6 (Wyithe 
& Loeb 2003a) . We use smoothed particle hydrodynamics 
(SPH) simulations to describe the cooling and dynamics 
of the gas for different choices of its total angular mo- 
mentum (or spin parameter), following the suggestion by 
Eisenstein & Loeb (1995) that low-spin systems would be 
more susceptible to the formation of a SMBH. If H2 cooling 
is suppressed inside these galaxies, then their gas will not 
cool below lO^K. When the temperature to which the gas 
cools is only somewhat lower than the virial temperature 
of the host galaxy (at which the pressure force balances 
gravity for the gas), we expect that fragmentation into 
small clumps will be avoided and the gas will tend to con- 
dense isothermally (at its temperature floor of ^ 10'*K) 
into large clumps. Such large clumps may then collapse 
to form a SMBH, possibly through an intermediate stage 
of a supermassive star. The viability of this scenario relies 
on the suppression of molecular H2 cooling, which when 
present is capable of cooling the gas to a temperature as 
low as 200 K. We therefore start our analysis in § 2 by 
considering the cosmic UV background expected from the 
first stars. Later on, in § 5 we carry out numerical sim- 
ulations to determine the level of UV radiation necessary 
for the photodissociation of H2 molecules during the early 
collapse phase of the above dwarf galaxies. 

The SPH code we use is the same as was used to sim- 
ulate the formation of the first massive stars (Bromm et 
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al. 1999, 2002), except that wc apply it to more massive 
galaxies where atomic; cooling is effective. Our goal is to 
find whether there is a direct route to the formation of 
SMBHs under these circumstances. Although our focus is 
on gas dynamics, we note that the first massive black holes 
(BHs) may have also formed out of the gravitational dy- 
namics of clusters of stellar-mass BHs (Larson 2000, 2002; 
Madau & Rees 2001; Schneider et al. 2002; Islam, Taylor, 
& Silk 2003; Volonteri, Haardt, & Madau 2003). 

Throughout this paper, wc assume a standard ACDM 
cosmology, with a total density parameter in matter of 
^lm = 1 — f^A = 0.3, and in baryons of JIb = 0.045. The 
Hubble constant is ft, = i?o/(100 km s"^ Mpc-^)=0.7, and 
the present-day power-spectrum amplitude is as = 0.9 in 
spheres of radius 8/i~^Mpc. 

2. THE COSMOLOGICAL CONTEXT 

Wc first outline the basic physical reason for why small 
protogalactic systems of virial temperature ^ 10*K and 
mass ~ IO^Mq, collapsing before the epoch of reioniza- 
tion at redshifts z ^ 10, provide intriguing sites for the 
formation of the first SMBHs. 

The direct collapse of a primordial gas cloud into a cen- 
tral compact object is made difficult by fragmentation and 
consequent star formation (Loeb & Rasio 1994). The abil- 
ity of the gas to fragment depends on the presence of an ef- 
ficient cooling mechanism. By a redshift of z ^ 10, the in- 
tergalactic medium (IGM) is expected to be enriched with 
heavy elements from the first generation of stars to a level 
of Z ^ 10-^-^ Zq (e.g., Gnedin & Ostriker 1997; Mackey, 
Bromm, & Hernquist 2003). Recent numerical investiga- 
tions have shown that such a metallicity is sufficient to en- 
able the gas to efficiently cool, fragment, and subsequently 
form stars (Omukai 2000; Bromm et al. 2001a; Bromm & 
Clarke 2002). It is plausible that in some regions of the 
IGM, the gas that collapses into a dwarf galaxy has not 
been enriched with metals in excess of ZcHt ~ 10~^'^Zq 
(e.g., Thacker, Scannapieco, & Davis 2002; Furlanetto & 
Loeb 2003). In this paper, we investigate the collapse of a 
virtually metal-free system of this type. 

It is well known that molecular hydrogen can effectively 
cool metal-free gas, and thus enable the formation of the 
first stars at 2: ~ 20 — 30 (e.g., Haiman, Thoul, & Loeb 
1996; Tegmark et al. 1997). Even in the absence of met- 
als, star formation is therefore expected to occur, and the 
subsequent negative feedback due to supernova explosions 
would prevent the assembly of large quantities of gas in 
the center of the shallow dark matter (DM) potential well 
that characterizes the first dwarf galaxies (e.g., Dekel & 
Silk 1986). Molecular hydrogen, however, is fragile and 
can readily be destroyed by photons in the Lyman- Werner 
(LW) bands, within the energy range 11.2 — 13.6 eV, via 
the two-step Solomon process (Stecher & Williams 1967) 

H2 + 7 ^ 2H. 

The intermediate stage involves an excited electronic state, 
H2 , from which a fraction of the subsequent decays end in 

the vibrational continuum of the ground state, resulting 
in the dissociation of the molecule. 



The question then arises whether H2 cooling can indeed 
be suppresscid in the pre-reionization dwarf galaxies con- 
sidered in this paper. These systems are rather massive 
compared to the ~ 10^ — lO^M© halos that host the for- 
mation of the first stars at z ^ 20. The gas might then be 
able to self-shield against the photo-dissociating LW back- 
ground (e.g., Glover & Brand 2001; Machacek, Bryan, & 
Abel 2001). Close to the epoch of reionization, however, 
a significant flirx in the LW bands {hi/ < 13.6 eV) is ex- 
pected. To estimate the LW flux, we first consider the flux 
of ionizing radiation just above the Lyman limit, 

+ he N^XpBjz = 10) 

Here, pb{z = 10) is the baryonic density at z ~ 10, 
X = 0.76 the mass fraction in hydrogen, and ma the mass 
of a hydrogen atom. We assume that A^-y 10 ionizing 
photons per baryon are required to reionize the universe 
(Wyithe & Loeb 2003a). Evaluating this expression for the 
cosmological parameters given above, we find: J21 ~ 40. 
Here, and in the remainder of the paper, UV fiuxes are 
normalized as J,y = J21 x 10"^"'^ erg cm~^ Hz~^ sr~^. 
The fiux in the LW bands just below the Lyman limit, 
could, however, be much larger than J+. Assum- 
ing that only a fraction, /esc, of the ionizing photons 
can escape from the star forming halos, we have: J21 — 
Jtil Icsc 4 X 10'^(/csc/0.01), where J^sc is expected to be 
low at high redshifts (Wood & Loeb 2000, and references 
therein)^. 

The strong UV background fiux estimated above makes 
it possible for H2 formation to be effectively suppressed 
close to the redshift of reionization. In the absence of 
molecular hydrogen, however, cooling can still proceed via 
atomic transitions in halos of mass (see Barkana & Loeb 
2001) 

/I . 2\ -3/2 

M^10«MQ(^^j • (2) 

The virial temperature, T^-^ ~ lO^K, in these more mas- 
sive halos allows for the very efficient cooling of the gas via 
lines of atomic hydrogen. Notice that in this case where 
the gas temperature is close to the virial temperature, the 
gas cloud as a whole can undergo collapse but it will not be 
able to fragment until high enough densities are reached 
so that the Jeans mass has declined sufficiently. A halo of 
total mass --^ lO^M© and a collapse redshift 2; ~ 10 cor- 
responds to a 2(7 peak in the random field of primordial 
density fiuctuations. 

In § 4, we describe mimerical simulations of such a 2 
cr peak. We consider the evolution with and without the 
presence of H2. To ascertain the importance of H2, we 
consider the limiting cases of a halo made of purely atomic 
hydrogen, and of a halo where H2 is allowed to form with- 
out any negative feedback. Subsequently, in § 5 we discuss 
simulations where an external LW background is included. 
We in particular assess the critical LW flux required to pre- 
vent molecules from forming. This required background 
flux may then be compared to the available flux at z ^ 10, 
as estimated in this section. 



strong background of photo-dissociating photons could also result from star formation in the dwarf galaxy itself (e.g., Omukai & Nishi 
1999; Oh & Haiman 2002). However, this would also result in the production of metals, probably leading to the rapid enrichment of the galaxy 
beyond .Zcrit- 
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3. NUMERICAL METHODOLOGY 

3.1. Large-scale Simulations 

The evolution of the dark matter and gas components 
is calculated with a version of TREESPH (Hernquist & 
Katz 1989), combining the SPH method with a hierar- 
(iiical (tree) gravity solver (sec Bromm ct al. 2002 for 
further details). Here, we briefly describe the additions 
to the code which are necessary for the investigation of 
zero-metallicity gas. These include a method to treat the 
radiative cooling of the gas together with a multi-species 
chemical reaction network, and a technique to create sink 
particles. The thermal evolution of the gas is governed by 
the equation: 

Du _ P Dp r - A 

where D/Di is the Lagrangian time derivative, P and p are 
the gas pressure and density, u is the specific internal en- 
ergy (in erg g~^), and T and A are the contributions from 
radiative heating and cooling, respectively (in erg cm^'' 
s^^). Since we are concerned with a system that formed 
prior to the epoch of reionization, we assume that there 
is no heating due to a photoionizing background (F = 0). 
Radiative cooling is due to lines of atomic hydrogen at gas 
temperatures T J> IC* K, and lines of molecular hydro- 
gen at lower temperatures. We have implemented the H2 
cooling function given by GaUi & Palla (1998). Finally, 
the first term on the right-hand side describes adiabatic 
cooling due to expansion or heating due to compression. 

Since radiative cooling to temperatures below that of the 
cosmic microwave background (CMB), Tcmb = 2.7 K(l -|- 
z), is thermodynamically not possible, we write the cooling 
term as 

A = A(T) - A(rcMB) . (4) 

For T < Tcmb, radiative cooling consequently turns into 
heating. This approximate treatment ensures that T > 
Tcmb, unless cooling proceeds via adiabatic expansion. 

The chemical reaction network comprises the 9 species 
H, H+, H-, H2, H^, e". He, He+, and He++, including 
the reactions given in Haiman et al. (1996). Our implicit, 
backwards-differencing method to solve the coupled set of 
rate equations is fast and accurate (see Bromm et al. 2002 
for test calculations). 

We have devised an algorithm to merge SPH particles 
in high-density regions in order to overcome the other- 
wise prohibitive time-step limitation, as enforced by the 
Courant stability criterion. To follow the simulation for a 
few dynamical times, we allow SPH particles to merge into 
more massive ones, provided they exceed a predetermined 
density threshold. More details of the merging algorithm 
are given in Bromm et al. (2002). 

3.2. Small-scale Simulations 

The sink-particle technique allows one to study the col- 
lective dynamics of multiple centers of condensation, such 
as merging between clumps and ongoing accretion over 
many dynamical timescales (e.g.. Bate, Bonnell, & Price 
1995; Bate, Bonnell, & Bromm 2003). It is important, 
however, to also be able to follow the collapse to in- 
creasingly high density. To accomplish this, we have car- 
ried out a complementary simulation in which we do not 



create sink particles. Instead, we focus on the highest- 
density region in one of the large-scale simulations, re- 
sample the density field with an increased number of SPH 
particles, and adopt a rapidly decreasing timestep accord- 
ing to Aisys cx l/^/Gpinax- Hcrc, /9,nax IS thc maximum gas 
density at a given instant, and Atgys, the system timestep, 
is the maxinium allowed time by which a particle is ad- 
vanced within the multiple-timestep scheme employed in 
the simulations (see Hernquist & Katz 1989). In doing so, 
the overall system does hardly evolve at all, while the run- 
away collapse of one small region proceeds on a sufficiently 
short timescale. The mass resolution of a simulation is ap- 
proximately 

«'" = (^)"» • ® 

where iVnoigh — 32 is the number of particles within a given 
SPH smoothing kernel, A^sph the total number of SPH 
particles, and Mb the total baryonic mass. To avoid nu- 
merical fragmentation, the Jeans scale has to be resolved 
(Bate & Burkert 1997): Mres < Mj. After the onset of 
gravitational instability in the large-scale simulation, the 
rapid increase in density leads to a violation of this cri- 
terion. At this point, a sink particle is created in the 
simulation. 

In our refined simulation, we follow a different compu- 
tational strategy (see also Bromm 2000). Our starting 
configuration is the large-scale simulation, stopped at a 
moment briefly before the central region undergoes run- 
away collapse. We now apply the following resampling 
procedure to the fluid within this region. Every SPH par- 
ticle in the original, unrefined simulation acts as a parent 
particle denoted by p, and spawns A'^ref child particles de- 
noted by k, where Nj-^f = 50 in this paper. The child 
particles are distributed according to the SPH smoothing 
kernel W{fk ~ ^p^hp), by employing a standard Monte- 
Carlo comparison-rejection method (see, e.g., Press et al. 
1992). Here, hp is the smoothing length of the parent 
particle. The velocity, temperature, and fractional abun- 
dances of the 9 species included in the reaction network 
are directly inherited from the parent particle, Vk = Vp, 
Tk = Tp, and yk = Xp for all fc, respectively. Finally, each 
child particle is assigned a mass of nrik = rUp/Nref- This 
procedure conserves linear and angular momentum well. 
Energy is also conserved well, although there arises a small 
artificial contribution to the gravitational potential energy 
due to the discreteness of the resampling. The resampling 
described in § 4.2 results in A^sph = 65,500 within the cen- 
tral high-resolution region, and the mass resolution is now 
Mrcs ^ 5OOM0, as compared to M^es ^ 20,000Mq in the 
original simulation. 

As we follow our simulation to densities in excess of 
riH ^ 10^ cm~^, opacity effects begin to modify the ther- 
mal behavior of the gas. For simplicity, we continue to use 
the optically thin rate for cooling due to the coUisional ex- 
citation of atomic hydrogen lines (Cen 1992). While this is 
clearly an oversimplification, we have verified that the re- 
sulting thermal evolution of the gas is very similar to what 
is found in a more sophisticated treatment of the cooling 
processes at high densities (Omukai 2001). 

Our technique of refining a coarser, parent simulation, 
and following the further collapse with increased resolu- 
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Table 1 

Parameters for the different Runs 



Run M 2vir A H2 LW nth 
(Mo) (cm-3) 

A 10^ 10 No No 10^ 

B 10^ 10 0.05 No No lO'' 

C 10* 10 Yes No 10* 

D 10^ 10 Yes Yes lO'' 



Note. — M is the total mass of the halo, Zvir the collapse redshift, A the spin pa- 
rameter, Titij the threshold density for the creation of a sink particle, H2 refers to 
the absence or presence of molecular hydrogen cooling, and LW denotes whether a 
photodissociating Lyman- Werner background is included or not. 



tion, is conceptually similar to the Adaptive Mesh Refine- 
ment (AMR) method which was originally developed by 
Berger & Oliger (1984). The AMR method has aheady 
been successfully applied to astrophysical problems (e.g., 
Truelove et al. 1998; Norman & Bryan 1999), but our ap- 
proach is one of the first attempts of implementing such 
a scheme within SPH (see also Kitsionas & Whitworth 
2002). 

3.3. External UV background 

To investigate the effect of an external UV background 
on the evolution of the gas, we have included the relevant 
photo-reactions into our chemical network. The photo- 
rates (in units of s~^) are given by 

^ = 4.^^^^, (6) 

where hvth is the threshold energy for a specific reaction. 
The cross-sections are given in Haiman et al. (1996). Fol- 
lowing Omukai (2001), we consider both a power-law spec- 
trum, oc appropriate for quasar sources, and a 
thermal Planck spectrum, Jj. cx B,^{T^), appropriate for 
stellar sources. In the latter case, we evaluate the photo- 
rates for two different values of the radiation tempera- 
ture: T, ~ 10^ K and 10^ K, corresponding to 'normal' 
Population II stars and to very massive Population III 
stars, respectively (Bromm, Kudritzki, & Loeb 2001b). 
For each spectrum, we normalize the flux at the Lyman 
limit with a normalization constant J21 , and assume zero 
flux at hv > 13.6 eV due to intergalactic HI absorption. 

To take into account the effect of self- shielding, we write 
the H2 photodissociation rate as /sdiss fx J21, /shield, with a 
proportionality constant that depends on the chosen spec- 
trum. For the shielding factor, we use the approximate 
expression /shield — min[l, (A^Ha/lO^^cm^^)^'^ '''''] (Draine 
& Bertoldi 1996). This formula is accurate only for a static 
mcditim, and it will grossly overestimate the H2 line opac- 
ity in the presence of large-scale velocity gradients. In our 
simulations, we encounter moderately strong bulk flows, 
with velocities of order the sound speed. We therefore con- 
servatively overestimate the effect of self-shielding. The H2 
column density is estimated from local quantities only as 
follows: N112 — O.l/Hj'^Hi'char- We define a local char- 
acteristic length, Lchar, such that the total baryonic mass 



would be contained in a sphere of uniform density nn, and 
radius Lchar- Thus, we have 

char 1^ 47rmH'^H / ^ ^ 

By comparing with the column densities that are actually 
realized in the simulations, we have verified the validity of 
this prescription. 

We neglect any contributions to the heating of the gas 
due to photo-reactions, since these are always negligible in 
the system studied here. We also neglect, for simplicity, 
the sawtooth-modulation of the background spectrum due 
to the IGM line-opacity in the Lyman series (Haiman et 
al. 1997). 

3.4. Initial Conditions 

Within the hierarchical ACDM model, the first luminous 
objects are expected to form out of high a peaks in the 
random field of primordial density fluctuations. The early 
(linear) evolution of such a peak, assumed to be an iso- 
lated and roughly spherical overdensity, can be described 
by the top- hat model (e.g., Padmanabhan 1993). We use 
the top-hat approximation to specify the initial DM con- 
flguration. In this paper, we investigate the fate of a 2 a 
peak of total mass IO^Mq, corresponding to 1.5 x lO^M© 
in baryons. On this mass scale, one can estimate the red- 
shift of collapse (or virialization) to be Zyir ~ 10. 

Our simulation is initialized at Zj = 100, by perform- 
ing the following steps. The coUisionless DM particles are 
placed on a cubical Cartesian grid, and are then perturbed 
according to a given power spectrum P{k) = Ak^, by 
applying the Zeldovich (1970) approximation which also 
allows to self-consistently assign initial velocities. The 
power-law index is set to n = —2.5 which approximately 
describes the spectral behavior on the scale ^ lO^M©. To 
fix the amplitude A, we specify the initial variance of the 
fluctuation amplitude 

C7f=^5]])fc" . (8) 

The summation is over all contributing modes, where the 
minimum wavenumber is given by the overall size of the 
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Cartesian box, and the maximum wavenumber by the 
Nyquist frequency. Choosing af ~ 0.01, the rms amph- 
tude of fluctuations at the moment of collapse is approxi- 
mately 



a{z = 10) 



1 + Zz 

1 + z 



1 



(9) 



This choice ensures that the substructure develops on a 
similar timescale as the overall collapse of the background 
medium. Next, particles within a (proper) radius of Ri = 
776 pc are selected for the simulation. The resulting num- 
ber of DM particles is TVdm = 17074. Finally, the par- 
ticles are set into rigid rotation and are endowed with a 
uniform Hubble expansion (see also Katz 1991). Angu- 
lar momentum is added by assuming a spin parameter 
A = L\E\^/'^/{GM^/^) = and 0.05, where L, E, and 
M are the total angular momentum, energy, and mass, re- 
spectively. The spin parameter is a measure of the degree 
of rotational support, such that the ratio of centrifugal 
to gravitational acceleration is given by ~ at virializa- 
tion. The; second value corresponds to the average spin 
parameter found in cosmological simulations (e.g., Barnes 
& Efstathiou 1987; Jang-Condell & Hernquist 2001). 

The collisional SPH particles (A^sph = 32768) are ran- 
domly placed to approximate a uniform initial density. 
The SPH particles follow the same Hubble expansion and 
rigid rotation as the DM particles. For the initial gas tem- 
perature at Zi = 100 we adopt the value of 200 K (Tegmark 
et al. 1997). The fractional free-electron abundance is ini- 
tialized as .Te = 4.6 X 10^^, and the hydrogen molecule 
abundance, in the runs where H2 formation is allowed, as 
= 2 X 10-*^ (Anninos & Norman 1996). 

In Table 1, we summarize the parameters of the large- 
scale simulations. We will discuss the initial state of the 
refined, small-scale simulation in § 4.2. We empasize that 
the initial particle setup, and in particular the realization 
of the DM fluctuations, is identical for all runs which only 
differ in the choice of spin, in the presence or absence of 
a soft UV background, and in whether H2 cooling is al- 
lowed or not. The runs with no H2 cooling correspond to 
the limiting case in which all the H2 has been radiatively 
destroyed by a soft UV background. In § 2, we have dis- 
cussed the possible emergence of such an efficient feedback 
following the build-up of the cosmic UV background close 
to the epoch of reionization. 

4. THE SIMULATIONS 

In this section, we only consider simulations where no 
external UV background is included, and defer the discus- 
sion of its effect to § 5. Wo carry out these simulations in 
two separate steps. First, we simulate the large-scale evo- 
lution on the scale of the host halo, ~ 1 kpc, ending with 
the formation of high-density clumps on a scale of ^ 1 pc. 
Once the mass resolution of the simulation becomes larger 
than the local Jeans mass, Mj <^ M^es, a sink particle is 
created. Consequently, the internal dynamics of a clump 
cannot be studied any further in these coarse-grain sim- 
ulations. Second, we follow the evolution of a clump to 
higher densities and smaller spatial scales, by refining the 
resolution in the vicinity of the clump to again resolve the 
Jeans mass. We begin with a discussion of the large-scale 
evolution. 



4.1. Large-scale evolution 

We first describe the run with no spin and with sup- 
pressed H2 formation (Run A) in more detail, and sub- 
sequently explore what happens under variations of these 
assumptions. In Figure 1, we show the initial configuration 
at Zi = 100, which is identical for all runs. The overdense 
region initially expands with the Hubble flow at a reduced 
rate, subsequently turns around at Zta, ^ 16, and eventu- 
ally collapses. 
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Fig. 1. — Run A: Initial configuration for top- hat collapse at 
Zi = 100. The halo has a total mass of 10*Mq, and is endowed 
with a Hubble expansion such that virialization occurs at Zvir ~ 10. 
Top row: The DM particles are perturbed from a regular grid ac- 
cording to P{k) oc Bottom row: The gas particles are placed 
at random, and comprise a mass fraction of 15%. Left panels: Face- 
on view. Right panels: Edge-on view. 



Dork Inciter 





-2000-lDOO 1000 2000 
X Position [pc] 



-2000-1000 1000 2000 
X Position [pc] 





-2000-1000 1000 2000 
X Position [pc] 



-2000-1000 1000 2000 
X Position [pc] 



Fig. 2. — Run A: Morphology at z ~ 13. The manner of presen- 
tation is the same as in Figure 1, but with a box size of ~ 5 kpc. 
The DM has developed significant substructure, and the baryons 
are just beginning to fall into the corresponding potential wells. 
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Fig. 3. — Run A: Morphology at 2 = 10.3. The convention in 
Fig. 1 is adopted for the rows and columns. The box size is r^ 2 
kpc. The DM is in the process of undergoing violent relaxation with 
the concurrent smoothing out of substructure. The gas has settled 
into the center of the DM potential well. 

Figure 2 shows the system at z ~ 13, briefly after 
turnaround. The DM component has developed a marked 
substructure in response to the imprinted perturbations, 
and the baryons have begun to fall into the deepest DM 
potential wells. Eventually, close to the redshift of virial- 
ization, Zvir — 10, the DM is undergoing violent relaxation 
(Lynden-Bell 1967), resulting in an approximate balance 
between kinetic and gravitational potential energy. 



X Position [pc] 

Fig. 4. — Run A: Central gas density at z = 10.3, in the case 
with zero initial spin. The box size is 200 pc. One compact object 
has formed in the center with a mass of 2.7 X IO^Mq and a radius 
of <1 pc. 

The baryons, on the other hand, dissipatively settle into 
the center of the DM halo (see Figure 3). Focusing on the 
innermost ^ 200 pc of the simulation box. Figure 4 dis- 
plays the central distribution of the gas density. At this 
stage in the simulation, a high-density sink particle of mass 



~ 3 X IO^Mq and radius 1 pc has been created. Sub- 
sequently, the sink particle grows in mass by continuous 
accretion from the surrounding diffuse gas. 
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Fig. 5. — Run A: Gas properties at z = 10.3. (a) Free electron 
abundance vs. hydrogen number density (in cm""^). (b) Hydro- 
gen molecule abundance vs. number density. In this run, we do 
not allow for the formation of H2, but compare with the corre- 
sponding panel in Fig. 7. (c) Gas temperature vs. number den- 
sity. At densities below ~ 1 cm" ^, the gas temperature rises be- 
cause of adiabatic compression until it reaches the virial value of 
Tvir ~ 10, 000 K. At higher densities, cooling due to atomic H keeps 
the gas nearly isothermal, (d) Jeans mass (in Mq) vs. number 
density. The Jeans mass reaches the resolution limit of the simula- 
tion, Mros — 20, OOOMq, for hydrogen number densities close to the 
merging threshold of rith = 10^ cm~'^. 

To gain further insight into the physics of the simulation, 
we consider the thermal and chemical properties of the gas, 
as illustrated in Figure 5. We show the free-electron abun- 
dance, temperature, and Jeans mass as a function of gas 
density for every SPH particle. This manner of presen- 
tation contains an additional dimension of information, 
concerning the overall timescale of evolution: when the 
evolution proceeds slowly, particles tend to pile-up in the 
respective plots, whereas only a few particles populate re- 
gions in parameter space where the evolution is fast. As 
can be seen in panel (d), a gas density of n ~ lO^cm"'^ 
and a Jeans mass of Mj ~ IO^Mq approximately mark 
the transition between a phase of slow and fast evolution. 
Physically, this transition corresponds to the onset of run- 
away collapse when the clump mass is Af ^ Mj. 

In Figure 6, we show the central gas configuration for the 
simulation where the initial spin of the system is nonzero 
(Run B). In this binary system of clumps has 

formed with a separation of ~ 10 pc. These two clumps in- 
dividually have properties similar to the one clump formed 
in Run A. The binary is rapidly drawn into the cen- 
ter of the halo by dynamical friction which operates on 
a timescale t^i ^ O.lr^Vc/iGMci) ~ 10^ yr (Binney & 
Tremaine 1987). Here, Vc is the circular velocity at radius 
r, and AIc\ is the clump mass. This timescale is short com- 
pared to the Hubble time at z 10, tn — 10* yr. Such 
a system of two compact objects is expected to efficiently 
radiate gravitational waves that could be detected with 
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the planned Laser Interferometer Space Antenna^ (LISA) 
(Wyithe & Loeb 2003b). 

Gas 







m 













-50 50 

X Position [pc] 



Fig. 6. — Run B: Central gas density at 2 ~ 10, in the case 
with an initial spin parameter of A = 0.05. The box size is 200 pc. 
Here, two compact objects have formed in the center with masses 
of 2.2 X 1O^M0 and 3.1 X 10^ M©, respectively, and radii <1 pc. 

Massive clumps, either in isolation or as a binary, are 
able to form in both runs A and B, which correspond to 
different values of initial spin. The result that such a com- 
pact object can form even in a halo with average initial 
spin (run B), is seemingly at odds with the earlier inves- 
tigation by Eisenstein & Loeb (1995), who have argued 
that only the lowest-spin cosmological perturbations can 
harbor seed BHs. The gas that is eventually incorporated 
into the two massive clumps in run B must have been able 
to efficiently lose angular momentum. It is crucial that 
the evolution in the A = 0.05 case leads to the formation 
of a binary system. Tidal torques can then transfer much 
of the angular momentum from the gas around each form- 
ing clump to the orbital motion of the system, similar to 
the case of present-day star formation in a clustered envi- 
ronment (e.g., Larson 2002; Bate et al. 2003). We plan 
to more fully address the complex physics of angular mo- 
mentum transport in virializing DM halos in future work. 

This efficient mechanism of assembling a large amount 
of gas into a compact (^ 1 pc) configuration crucially de- 
pends on the suppresion of star formation in the collapsing 
halo. In Run C, which otherwise has the same initial con- 
ditions as Run A, we now allow for the formation of hy- 
drogen molecules. Figure 7 (which should be compared to 
Fig. 5 of Run A), shows that the thermal evolution of the 
gas in this case proceeds very differently from Runs A and 
B. Here, the gas is able to cool to temperatures T ^ 200 K 
due to the presence of II2 . Provided that II2 formation is 
not suppressed, the properties of primordial, metal-free 
gas are rather similar in systems of mass ~ 10^ Mq, and 
in the smaller halos of mass IO^Mq that are predicted 
to host the very first stars at z ~ 20 — 30 (see Fig. 10 in 
Bromm et al. 2002). 

The thermal evolution is clearly reflected in the mor- 
phology. Figure 8 shows the projected gas density. Due 

®See http://lisa.jpl.nasa.gov/ 



to efficient cooling, the gas can readily fragment and sub- 
sequently undergo runaway collapse to form stars already 
shortly after turnaround. This star formation activity will 
tend to prevent the assembly of gas in compact clumps in 
two ways. First, the Population III stars are expected to 
be short-lived, with lifetimes of ~ 3 x 10^ yr, and they 
will exert a strong negative feedback on their surround- 
ings upon their death as supernovae (SNe). Second, the 
SNe will disperse the nucleosynthetic products from the 
first generation of stars into the remaining gas of the halo. 
Henceforth, the enriched gas will be able to cool and form 
stars regardless of whether H2 is present or not. Since 
the photodissociation of H2 is now no longer able to limit 
the gas condensation (see § 2), the star formation rate is 
expected to increase substantially (e.g., Nishi & Tashiro 
2000; Bromm & Clarke 2002). 




-2Q246S -2 02488 

log Hh [crn~^] loq n» [crrt""*] 



Fig. 7. — Run C: Gas properties at z ~ 12. (a) Free electron 
abundance vs. hydrogen number density (in cm~'^). At densities 
exceeding njj ~ 10^ cm"'', recombination is very efficient, and the 
gas becomes almost neutral, (b) Hydrogen molecule abundance vs. 
number density. In this run, we do allow for the formation of H2. 
The thermal behavior of the gas is markedly different from Runs A 
and B, where molecule formation was suppressed, (c) Gas temper- 
ature vs. number density. In contrast to Runs A and B, cooling due 
to H2 drives the temperature down at higher density, to values of ~ 
200 K. (d) Jeans mass (in Mq) vs. number density. The Jeans mass 
reaches the resolution limit of the simulation, Mres — 20, OOOMq , 
for densities close to the merging threshold of = 10* cm""'. 

The massive clumps formed in runs A and B are evi- 
dently not yet BHs. To examine the subsequent dynamics 
of the clump, we have resimulated the evolution of a typi- 
cal clump with higher resolution in order to reach smaller 
spatial scales. We report on this simulation next. 

4.2. Small-scale evolution 

As a representative example, we have selected Run A for 
initializing the refined simulation. We focus on the central 
region, as shown in Figure 4, and resample the gas within 
a radius of ~ 5 pc from the density maximum. The ini- 
tial moment of the resampled run is chosen to be somewhat 
earlier than in Figure 4, corresponding to z ~ 10.6, slightly 
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before a sink particle was created in the coarse-grain sim- 
ulation. In the resampled region, the baryons comprise a 
mass of ~ 6 X IO^Mq, and are now the dominant compo- 
nent, compared to only IO^Mq in dark matter. 

After ^ 10'' yr, we reach the end-state of our refined 
simulation, at which point the Jeans mass again becomes 
comparable to the resolution limit. The central gas cloud 
is in a state of free-fall, and we do not see any signs of fur- 
ther sub-fragmentation. We have verified that fragmenta- 
tion is not artificially suppressed in the resampling process. 
To this extent, we have carried out a fiducial simulation 
that does lead to fragmentation at some point in the evo- 
lution, and a comparison calculation where we resample 
the fluid before fragmentation occurs. We find that the 
resulting fragmentation pattern is very similar in the two 
cases, as is physically expected. 

The innermost region, of size ^0.1 pc, comprises ^ 
10'' Mq in gas with densities in excess of 10^ cm^'^. If 
we were to continue the simulation further in time, the 
amount of gas residing inside the central ^0.1 pc would 
rapidly increase. Indeed, the mass of the sink particle 
formed in the large-scale simulation of Run A, of order a 
few IO^Mq , is indicative of the total baryonic mass that 
will end up in the central compact object. At this stage, 
this object is characterized by a ratio -Erot/I-Bgravl — 0.5, 
where -Eiot and -Egiav are the rotational and gravitational 
energies of the central compact object, respectively. The 
angular momentum present in the central clump must have 
arisen through torques from the clumpy DM distribution 
during the relaxation process, as run A has zero initial 
spin. Again, we plan to revisit this issue in our planned 
work on angular momentum transport in collapsing DM 
halos. 




Fig. 8. — Run C: Morphology at z ~ 12. Shown is the projected 
gas density in the x-y plane. The box size is ~ 3 kpc. Here, the gas 
readily fragments and subsequently undergoes runaway collapse to 
form stars due to the efficient cooling provided by H2. 

The crucial question now is: What is the further fate of 
the central object? In particular, do we expect the cloud 
to fragment, and eventually form a central stellar cluster? 
Or else, will fragmentation be inhibited, as in the earlier, 
large-scale evolution described in the previous section? In 



the latter case, the additional question arises (Loeb & Ra- 
sio 1994) whether the gas cloud will continue to collapse 
directly to a black hole or else settle temporarily into a 
pressure-supported configuration of a supermassive star 
(see §17 in Shapiro & Tcukolsky 1983)? We can gain im- 
portant physical insight into these questions by comparing 
three important timescales in the problem, namely: the 
free-fall time, iff, during which dynamical equilibrium is 
established; the cooling time, icooi, during which the gas 
may radiate its thermal energy as it contracts and heats 
to higher temperatures; and the viscous time, fyis, dur- 
ing which angular momentum is transported. The cooling 
time remains relatively short and allows continued col- 
lapse until the cloud becomes optically thick to the ra- 
diation it produces. Spherically-symmetric calculations 
indicate that when the particle density in the contract- 
ing gas clump rises above ~ 10^^ cm~'^, the ionization 
fraction increases sharply by 7 orders of magnitude to a 
value Xe <; 0.1 (see Fig. 4 in Omukai 2001). A clump 
of mass ~ 2 X IO^Mq achieves this dense phase when its 
radius R shrinks to (i?/10^^ cm) = Riq <^ 0.5, about two 
orders of magnitude smaller than the value at the end 
of our simulations. When Riq <^ 1, the optical depth 
of the cloud to Thomson scattering by free electrons is 
r ~ lO^R^Q and so the cooling time is given by the pho- 
ton diffusion time icooi ~ tR/c ~ 10^ yr Riq - The free 
fall time, tg ~ {Gp)~°-^ ~ 1 yr ^le^, is much shorter than 
this cooling time and so the gas may settle into a radiation- 
pressure supported configuration resembling a rotating su- 
permassive star (SMS). The viscous time is shorter than 
the cooling time, tvis ^ 4 x lO" yr Rli?{Pg^/ Ptot)a-^ , 
where a ^ 1 is a dimensionless measure of the viscosity 
coefficient (Shakura & Sunyaev 1973; §14.5 in Shapiro & 
Teukolsky 1983) in terms of the total (gas-|-radiation) pres- 
sure, Ptot = (Pgas + -Prad) > -Pgas- The resulting viscous 
dissipation will likely heat the optically-thick gas close to 
its virial temperature since the rotational energy is sub- 
stantial at the end of our simulation. Combined with ad- 
ditional heating from adiabatic compression, the collapse 
may be halted as soon as the gas becomes highly optically- 
thick and is unable to cool efficiently. We therefore con- 
clude that if the simulated gas clouds continue to shrink 
during the early optically-thin (efficient cooling) regime by 
about two orders of magnitude in radius (as they do earlier 
in the simulations), then they will likely form a rotating 
SMS. 

Recent fuUy-relativistic calculations of the evolution of 
rotating SMSs have predicted the collapse into a massive 
BH (Baumgarte & Shapiro 1999; Shibata & Shapiro 2002). 
Under a wide range of initial conditions, a substantial frac- 
tion of the mass of the SMS (~ 90%), is expected to end 
up in the BH. 

5. EFFECT OF A PHOTODISSOCIATING BACKGROUND 

We now return to the question: What is the required 
level of the LW background to prevent H2 from forming 
when self- shielding is taken into account? To answer this 
question, we have carried out a series of simulations where 
we allow H2 to form self-consistently in the presence of a 
soft UV background. We have determined Jcrit, the crit- 
ical flux above which II2 is not able to form throughout 
the simulation volume, for the three spectra described in 
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§ 3.3. For both the power-law spectrum, as weU as the 
thermal spectrum with T^, = K, wc find: Jcin > 10^. 
This value is larger than the expected background close to 
reionization (see § 2), and it would be difficult to prevent 
H2 cooling for these two spectra. 

A much lower background level, however, suffices in the 
case of a thermal spectrum with = 10"* K, for which we 
find: Jcrit ^ 10^, close to the predicted UV background at 
2 ~ 10. In Figure 9 we show the timescales that are rele- 
vant to understand the abundance of H2 in the simulation 
with J21 = 10"^. It is evident that even when self-shielding 
is included, the H2 formation time, t[orm, is longer than 
the H2 destruction time, tdest throughout the simulation. 
Once the gas has reached a density of nn ^ 10^ cm"-^, 
collisional dissociation becomes very effective in destroy- 
ing H2 at a gas temperature of lO^K. The ever-increasing 
column densities, and the resulting level of self-shielding, 
have therefore no effect on the molecule abundance at den- 
sities high enough for collisions to become important. 

Using an idealized one-zone model to describe the dy- 
namics of a collapsing primordial cloud, but otherwise im- 
plementing the relevant radiative processes in a sophisti- 
cated way, Omukai (2001) has determined the critical UV 
flux in the case of a Planck spectrum with = 10** K 
to be: J21 ~ 10^, in good agreement with our numerical 
simulations. 




log n„ [cm 

Fig. 9. — Run D: Important timescales at z = 10.6 for a back- 
ground stellar spectrum w ith T, = 10* K. Solid line: Free-fall time 
vs. hydrogen number density (in cm~^). Dots: H2 formation and 
destruction times vs. tin for every SPH particle in the simulation. 
The LW flux is normalized to J21 = lO"' at the Lyman limit. The 
destruction time increases with density up to njj ~ lO'^ cm^'* due to 
the effect of self-shielding. H2 is never formed efficiently, however, as 
iform > tdost throughout the simulated volume. At rtn > 10'^ cm^'\ 
H2 is readily destroyed via collisions, and tjest decreases again with 
density. 

Both these estimates crucially depend on the presence of 
a stellar-like radiation background. The formation of the 
first quasars, therefore, might well require an earlier epoch 
of star formation. The old puzzle of whether quasars or 
stars were the first luminous objects in the universe (e.g.. 
Silk & Rees 1998) would then be answered in favor of the 
latter. 



6. SUMMARY AND CONCLUSIONS 

Our numerical simulations show that metal-free dwarf 
galaxies at redshifts z ^ 10 whose cooling is dominated by 
atomic transitions (with a virial temperature just above 
lO^K) tend to form a central condensation consisting of 
one or two clumps which contain ^ 10% of the total bary- 
onic mass of the galaxy. As long as H2 formation is sup- 
pressed, these massive clumps do not fragment but rather 
cool and continue to collapse isothermally at a temper- 
ature of 10*K. We have simulated the collapse to the 
point where ^ IO^Mq of gas condense to a scale ^ Ipc. 
At the end of our simulation the clump maintains a nearly 
free-fall collapse and shows no signs of fragmentation. We 
expect the collapse to continue until the gas cloud becomes 
optically thick to Thomson scattering and its cooling time 
is much longer than the infall time. At this stage the col- 
lapse may be halted due to viscous dissipation of rotational 
energy, and a rotating SMS is likely to form. After shed- 
ding a small fraction of its mass, the SMS will collapse to 
a SMBH (Baumgarte & Shapiro 1999; Shibata & Shapiro 
2002). 

The possible existence of a supermassive star as an 
intermediate stage in the formation of the first quasars 
can be tested by direct observations. The SMS spec- 
trum is expected to be close to thermal at an effective 
temperature of ^ 10"'' K (Bromm ct al. 2001b). The 
SMS mass can be inferred from its bolometric luminos- 
ity, given that the latter equals the Eddington value, 
LEdd = 1-4 X 10*^ erg s''^ {M/ 10^ Mq). The emission red- 
shift can be easily derived from the Gunn-Peterson trough 
produced by Lya absorption of the neutral IGM (Gunn 
& Peterson 1965). Our model predicts that the emission 
spectrum of the SMS would show only hydrogen and he- 
lium lines but no metal lines (such as the broad emis- 
sion lines observed for quasars at lower redshifts), since 
metal enrichment would signal small scale fragmentation 
into stars and would in turn enable the production of 
molecules on dust grains. However, all the above char- 
acteristics are common to clusters of Population III stars 
(Bromm et al. 2001b). To distinguish between a SMS and 
the extended image of a star cluster requires high angu- 
lar resolution, possibly achievable with the James Webb 
Space Telescope^. Alternatively, it is possible to make use 
of the high probability of gravitational lensing by interven- 
ing galaxies out to z ~ 10 (Barkana & Loeb 2000; Wyithe 
& Loeb 2002a,b; Comerford, Haiman, & Schaye 2002). In 
particular, gravitational microlensing by stars within the 
lens galaxies (Wyithe & Loeb 2002b) is only possible for 
a SMS and not for a star cluster whose extent is much 
larger than the Einstein radius of a solar-mass star at a 
cosmological distance, ~ 10~^ pc. 

Our calculations show that fragmentation into stars is 
suppressed inside the above dwarf galaxies only if II2 cool- 
ing is negligible. The photodissociation of H2 molecules re- 
quires a background flux of UV photons below 13.6eV (to 
which the neutral IGM is nearly transparent) that could 
possibly be produced by stars at ^; ^ 10, prior to reioniza- 
tion (see §2). 

What is the comoving density of the first generation of 
massive BHs? At z ~ 10, the comoving density of halos 



See http://ngst.gsfc.nasa.gov/ 
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with a virial temperature ^ IC^K is ^ 5 Mpc (Shcth & 
Tormen 1999). Requiring that the total mass density of 
massive BHs would not exceed the value observed in the lo- 
cal universe (Yu & Tremaine 2002), ~ 2.5 x IO^Mq Mpc-^, 
wc infer that the average BH mass per halo must be 
^ 5 X IO^Mq, much smaller than the characteristic mass of 
the gas clumps in our simulations. This follows naturally 
from our expectation that only a small fraction of all halos 
are both metal poor and exposed to a high level of UV flux 
during their entire history (see §2 and §5). Even though 
a substantial fraction of the IGM volume is not enriched 
with metals at 2; ^ 10 (Thacker et al. 2002; Furlanetto & 
Locb 2003), most 2cr halos may have been enriched by stars 
that were formed inside their hierarchical building blocks 
at an earlier cosmic time, when the UV background was 
much lower. The metals may then facilitate H2 formation 
on dust particles as well as efficient cooling through metal 
emission lines (both atomic and molecular), resulting in 
further fragmentation. In such a case, only a small mi- 
nority of all halos forms by direct accretion of pristine gas 
from the IGM, as we have assumed in our simulations. 

Alternatively, the final BH may contain a small frac- 
tion of the initial mass of its parent clmnp of gas. Since 
we did not follow the collapse of each clump down to its 
Schwarzschild radius, we cannot be certain of the final BH 
formation efficiency. It is possible that hydrodynamic or 
radiative feedback from a growing BH seed at the cen- 
ter limits the final mass that it obtains, so that eventu- 
ally most of the mass of its parent clump is expelled in a 
wind due to energy release near the center (e.g., Haehnelt, 
Natarajan, & Rees 1998; Silk & Rees 1998). The corre- 



lation between BH mass and the velocity dispersion of its 
host galaxy as observed in the local universe (Merritt & 
Ferrarese 2001; Tremaine et al. 2002) would imply an ex- 
ceedingly low efficiency (^ 10~*) for the BH formation 
process in the dwarf galaxies under consideration (Wyithe 
& Loeb 2002c). However, this phenomenological correla- 
tion is measured in metal-rich galaxies where molecular 
cooling and star formation are abundant, and it may not 
hold under the unusual physical conditions inside the first 
generation of metal-free galaxies where the seeds for the 
present-day population of BHs were produced. 

An early formation of SMBHs would lead to a con- 
siderable rate of detectable bursts of gravitational radia- 
tion from the coalescence of BH binaries at high redshifts 
(Wyithe & Loeb 2003b). In addition to the commonly 
discussed origin of binaries in galaxy mergers (a possi- 
ble example of which has recently been observed in the 
nearby galaxy NGC 6240 by Komossa et al. 2003), we 
have demonstrated that BHs may form in binaries to start 
with, just as stars do. 
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able to us a version of TREESPH. We also thank Rennan 
Barkana and the anonymous referee for valuable comments 
on the manuscript. This work has been supported in part 
by NSF grants AST-0071019, AST-0204514. VB thanks 
the Institute for Advanced Study for its hospitality dur- 
ing the completion of this work. AL acknowledges support 
from the Institute for Advanced Study at Princeton and 
the John Simon Guggenheim Memorial Fellowship. 
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